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RELATIVE ROLE OF STARS AND QUASARS IN COSMIC REIONIZATION 
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ABSTRACT 

We revisit the classical view that quasar contribution to the reionization of hydrogen is unimportant. When 
secondary ionization are taken into account, in many plausible scenarios for the formation and growth of super- 
massive black holes quasars contribute substantially or even dominantly at z > 8, although their contribution 
generally falls below that of star-forming galaxies by z = 6. Theoretical models that guide the design of the first 
generation of redshifted 21 cm experiments must, therefore, substantially account for the quasar contribution 
in order to be even qualitatively accurate. 

Subject headings: cosmology: theory - cosmology: large-scale structure of universe - galaxies: formation - 
galaxies: intergalactic medium 
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1. INTRODUCTION 

Modeling theoretically the properties of the sources that re- 
heated and re-ionized the intergalactic medium at the end of 
the cosmic Dark Ages is difficult, because, at present, our 
knowledge of the nature of ionizing sources is highly in- 
complete. The con ventional view, held over a decade since 
the original work of Shapiro & Giroux' ('1987"); Shapiro et al.' 
([1994); Giroux & Shapiro (1996a); Madau et al. (1999), is 
that the ionization at the highest redshifts is dominated by 
the UV radiation from star-forming galaxies, and that the 
non-thermal contribution from quasars builds up only later, 
at z < 4, to dominate the double ionization of helium (which 
requires photons of energies exceeding 54.4 eV). 

This conventional view, however, is incomplete. In partic- 
ular, it ignores an important physical process of "secondary 
ionizations" by energetic ionizing photons. It is a well-known 
physical fact that an ionizing photon with the energy E in 
excess of ^ lOOeV is capable of ionizing more than a sin- 
gle atom - the excess photon energy of E - 13.6 eV is de- 
posited in the electron, which is cap able of ionizing one or 
more additional atoms in its vicinity (IShull & van SteenbergI 
[T98llValdes & Ferrarall2008l) . 

These secondary ionizations can substantially increase 
the relative contribution of quasars to the reionization of 
the universe. For example, a 1 keV X-ray photon can 
ionize up to ^ 25 hydrogen atoms. While this important 
effect has been included in many of the prior work on 
hydrogen reionization dOiroux & Shapiro 1996b; Oh 2001; 
Venkatesan et al. 2001; Ricotti et al. 2002; Machacek et al. 
2003; Ricotti & Ostriker 2004; Diikstra et al. 2004; 
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hydrogen reionization in light of recent improvements in the 
values of cosmological parameters and new developments in 
our understanding of the growth of supermassive black holes 
at high redshifts. 

There are two possible approaches to modeling reioniza- 
tion. In the first one, a model for the emission, propagation, 
and absorption of ionizing radiation from various categories 
of sources is constructed. Ultimately, such a model must in- 
volve a cosmological numerical simulation that resolves Ly- 
man limit systems with sub-kpc resolution and clustering of 
galaxies and quasars on ^ 100 Mpc scales. Such a large uni- 
form dynamic range is not yet feasible in modern cosmolog- 
ical simulations, so this type of modeling inevitably suffers 
from the unknown systematic biases, but, at the same time, it 
offers the most complete model of cosmic reionization. 

A second, much simpler (and, therefore, much more lim- 
ited) approach consists in ignoring the absorption and propa- 
gation of ionizing photons, and restricting a theoretical model 
to only counting the number of ionizations per atom. A ma- 
jor (but not the sole) limitation of this approach is that the 
criterion for reionization - the required number of ionizations 
per atom - is unknown. While reasonable estimates for this 
number can be constructed (e.g. lGnedinl l2008. and discussion 
there), it, at the very least, must exceed unity^. 

In this brief paper, we adopt the second, simple approach, 
to couple theoretical expectations for both stellar sources and 
quasars, and to consider their relative contribution to the to- 
tal number of ionizations per hydrogen atom. In line with 
our limited goals, we restrict our attention only to two main 
astrophysical types of ionizing sources, although we remain 
acutely aware that other, more exotic sources - including the 
energetic photons and electrons from dark matter annihilation 
- can provide a substantial or e ven dominant contributi on to 
the total ionization budget (e.g. lBelikov & Hoopeiil20()9L and 
references therein). 

2. METHOD 

We evolve the population of two distinct types of seed 
massive black holes (MBHs): either "small seeds", de- 
rived from Po pulation III remnants (Madau & Rees 200il 
IVolonteri et alJl2003b , or "large seeds", de rived from gas- 
dynamical collapse in metal-free galaxies ( iBegelman et alj 

^ Or, even more precisely, about 0.97, as a few percent of gas remains 
neutral after reionization, locked in Lyman limit and IDamped Lya systems. 
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l2006tlLodato & NataraiaDl2006[) . We adopt one single forma- 
tion scenario for each realization of the Universe we consider, 
therefore each model contains only "small seeds" or "large 
seeds". 

In the Population III remnants model, seed MBHs form 
with masses Wseed ^ fewx 10^ M© (see Fig. 1, top left panel), 
in haloes collapsing at z > 15 from rare 3.5 -a peaks of the 
primordial density field jVolonteri et alj|2003h . 

In the "large seeds" scenario, massive seeds with M « 
10'*Mq can form at high redshift, when the intergalac- 
tic medium has not been significantly enriched by metals. 
Here w e refer to lBegelman et alJ (l2006l) : lLodato & NatarajanI 
(iOOi), for more details of the physical MBH formation 
model. Seeds form in gravitationally unstable pre-galactic 
disk with primordial composition, in halos with virial tem- 
perature ^ lO'* K, cooled mainly by atomic hydrogen. The 
stability of gaseous disks depends on two parameters, the halo 
spin parameter Aspin, and the fraction of baryonic matter that 
ends up in the disk, = (flM / ^ h){M disk /Mh). A maximum 
spin parameter Ajpin.max exists, for which a disk is unstable as 
a function of the fraction of baryons forming the disk, i.e., for 
every /j, disks are stable for Aspin > Aspin.max- We here assume 
(fd, Aspin,max)=(0.2,0.2). The mass of the forming MBH seeds 
is set by the joint characteri stics of the gas flow and of the evo- 
lution of the collapsed gas dBegeknan et alj|2007l) . The mass 
function of MBH seeds peaks at lO'^M© (see Fig. 1, top right 
panel). 

We study MBH evolution within dark matter halos via 
a Monte-Carlo algorithm based on the extended Press- 
Schechter formalism. The population of MBHs evolves along 
with their hosts according to a "merger driven scenario", as 
described in Volonteri et al. (2003; 2009). An accretion 
episode is assumed to occur as a consequence of every major 
merger (mass ratio larger than 1:10) event. During an accre- 
tion episode, each MBH accretes an amount of mass, AM = 
9 X 10^ Mo((T/200km/s)'*, that scales with the Mbh-ct* rela- 
tion of its hosts (see Volonteri & Natarajan 2009). Accretion 
starts after a dynamical timescale and lasts until the MBH has 
accreted AM. 

We model the accretion rate onto MBHs in two ways. As 
a baseline model, we as sume that accretion p roceeds at the 
Eddington rate (see also Salvaterra et alj|2007h . In a second 
case, we model the accretion rate during the active phase from 
the extrapolation of the empirical distribution of Eddingto n 
ratios, A = log(Lboi/iEdd), found in iMerloni & Hein"3 (|2008|) . 
We adopt a fitting function of the Eddington ratio distribution 
as a function of MBH mass and redshift (Merloni 2009). We 
are here extrapolating such a model at much higher redshifts 
and lower MBH masses than originally intended. We caution 
readers in taking the results of this model face value. The main 
goal of our exercise is to probe possible sensible ranges for the 
accretion rates on MBHs. 

At each timestep of our Monte Carlo simulation of halos 
merger trees, we calculate the average energy density emitted 
during each timestep by the MBH population as follows: 

AU={e)Ap,,,c\ (1) 

where Apacc is the total mass density (in comoving units) ac- 
creted by MBHs within the timestep, and (e) is the average ra- 
diative efficiency, which we assume depends solely on MBH 
spins for radiatively efficient accretors. We evolve MBH 
spins accor ding to two simple models: cohere nt accretion 
jVolonteri et al. 2005) and chaotic accretion ( King & Pringld 
120071) . These two models lead to rapidly spinning MBHs 
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Fig. 1. — Mass function of MBHs at different cosmic times. Top to 
bottom: 18 < z < 20; 16 < z < 18; 14 < z < 16; 12 < z < 14; 10 < z < 12; 
8 < z < 10; 6 < z < 8. Left panels: small seeds. Right panels: large seeds. 
Solid histogram: Eddington accretion; dot-dashed histogram: Merloni-Heinz 
accretion. 

("high spins", (e) = 0.2), and slow rotators ("low spins", 
(e) = 0.06) respectively. We assume that MBHs are born non- 
spinning. 

The energy density, U, depends on both the accretion rate 
(in Eddington units, fEdd = lO'^) and e. Assuming constant 
fEdd and e, at a given time t: 

U(t) oc eMseed {e^"" ^ ^ - 1 ) (2) 

where Mseed is the average mass of MBH seeds in each sce- 
nario, and T = 0.45 Gyr is the Eddington timescale. We can 
see that if f ^ r: U oc fnddi^ -e)- Within our evolutionary 
scheme e is determined for each MBH, while we have used 
a single, average, (e) in calculating U. Since it takes a few 
Myr for MBHs in the coherent accretion model to spin-up to 
large spins (i.e. for the radiative efficiency to increase from 
0.06 - Schwarzschild hole, to 0.2 - spin= 0.9) this scheme 
slightly overestimates the radiative output from MBHs in the 
'high spin' case at z > 11. 

We further assume that a fixed fraction /uv of the bolo- 
metric power radiated by high-redshift quasars is emitted as 
hydrogen-ionizing photons. The number of ionizing photons 
scales as fuy/E^, where is the mean photon energ y (see 
section 2). We refer the reader to Madau et al. (2004) for a 
thorough discussion about the quasar spectra; in the follow- 
ing we assume, conservatively, that /uv = 0.2. 

The total comoving energy density of ionizing radiation 
emitted by the growing MBH population is then: 

U^(z)= ^ fmm (3) 

i=OJ(z) 

where the sum is over the timesteps between our starting red- 
shift (Zmax = 20, timestep 0), and the timestep, j{z) correspond- 
ing to the redshift of interest, z- 

We have also implemented a more pessimistic case fo r 
MBH accretion, loosely based on lMilosavljevic et alJ (l2008l) . 
where the accretion rate has been fixed to 30% of the Ed- 
dington accretion rate. In this scenario, the yield in ionizing 
radiation is lower than for the case of 100% Eddington rate. 
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but ionization histories in this scenario fall within the range 
spanned by our other models; they, thus, do not change any of 
our conclusions and we do not consider them further in this 
paper 

Our models are consistent with the constraints from the soft 
X-ray background (Dijk stra et al. 20 04b). Our predicted pop- 
ulation of high redshift AGN would account for almost 5% of 
the measured (0.5-8 keV) bac kground, or ~ 25% of the un- 
resolved one jSalvaterra et al.1 2005, 2007). Models are also 
consiste nt with t he bolometric luminosity function of quasars 
at z > 4 dHgpkins et al. 2007). "Large seeds" are preferred as 
our fiducial model, based on a better agreement with the bolo- 
metric luminosity function at z ~ 4, however we note that in 
this paper we pr obe much hig her redshifts than those probed 
by the study of Hopkins et al. (2007), hence we consider the 
match with the luminosity function at lower redshifts as a 
weak constraint. 

For the stellar contribution to the total ionizing background 
we use the extrapolation of the o bserved UV luminosity func- 
tions of high redshifts galaxies jBouwens et alj|2007l l2008b 
and the assumed value for the relative (to the amount of es- 
caped UV light at 1000) escape fraction of ionizing radia- 
ti on. The c omple te details of the methodology are described 
in iGnediiil (l2008h . Here we only briefly repeat that the total 
mass-to-light ratios of higher redshift galaxies are computed 
in a given cosmology by matching the observed spatial abun- 
dance of galaxies of a given luminosity to the theoretically 
computed abundance of halos of a given mass, so that 

ni> Luv) = ni> Mtot) (4) 

(and, optionally, an additional factor that accounts for 
"bursty" star formation rate in high redshift galaxies can be 
incorporated in eq. ID; reasonable values for that factor make 
insignificant impact on our results). 

The obtained thus mass-to-light ratio can be extrapolated to 
higher redshifts. Since the mass-to-light ratio is a weak func- 
tion of redshift, the uncertainty of such an extrapolation does 
not dominate the final uncertainty of our estimate for the ion- 
izing emissivity from galaxies. Instead, the final uncertainty 
is dominated by the uncertainty (both observational and theo- 
retical) on the ado pted value of the relative escape fraction for 
ionizing radiation dOnedinllfOO S). 

The combined (stellar plus quasar) contribution to reioniza- 
tion can then be estimated as 

^'^/''=^'^/-^^^^«'^^^T4:Sv;?:' 

where N^/„ is the number of ionizations per atom, «« = (!- 
0.75Yp)ni, « 2.0 x 10"^ cm"-' is the comoving number den- 
sity of atoms of hydrogen or helium (we assume that most 
of helium is only singly ionized during hydrogen ionization), 
and 14.4 eV is the mean ionization energy per atom. The 
first term in this equation accounts for the contribution from 
stars, the second one includes primary ionizations by ioniz- 
ing photons from quasars with the mean photon energy E^, 
and the last term accounts for secondary ionizations from en- 
ergetic ionizing photons. For the gas uniformly ionized to the 
ionization fraction x, the fraction of radiation energy density 
/sT going into secondary ionizations has been computed by 
IShuU & van Steenberd (11985) . and their results c an be conve- 
niently fitted by a simple but accurate formula dRicotti et al.l 



l2002h 

/si(x)«0.35(l-x«-4)'-«-1.77 (^^y'xO-2(l-x«Y. 

(6) 

In reality, of course, the universe is not uniformly ionized, so 
the quantity x in equation (|5j is an effective value Xeff, such 
that 

and the average is mass-weighted. Thus, Xeff is not a mass- or 
volume-weighted average of the cosmic ionization fraction. 
However, jCeff = for the fully neutral and Xetf = 1 for the fully 
ionized universe. 

The quantity N^/a.* has been computed in IGnediiil ( l2008h 
and is not discussed here. The mean energy of ionizing pho- 
tons from quasars, E~^, remains a parameter of our model. 
The exact value of this parameter obviously depends on 
the spectral shape of the quasar energy distribution. If we 
assume a classic multicolor disk spectrum up to kT^^y, 
1 keV(MBH/ M0)-'/4 (Shakura & Sunyaev 1973), and a non- 
thermal power-law component with spectral slope cx t^"", 
with a ~ 1 at higher energies, we find ~ 300 eV for Mbh = 
10^ Mq. The spectrum is harder/softer for smaller/larger 
MBH masses. As a fiducial value, we adopt E~^ = 300 eV (cfr. 
Fig 1), and we investigate the sensitivity of our results to the 
precise value of this parameter in the next section. 

Equation ^ cannot be evaluated without an equation for 
jCeff as a function of time. A complete computation of Xeff(r) 
requires a sophisticated three-dimensional modeling of the 
transfer of ionizing radiation throughout the inhomogeneous 
gas density distribution in the universe. Since such modeling 
is well beyond the scope of this paper, we instead introduce 
the following ansatz forxeff(r): 

jceff = min( 1 , fn^xN^la), (7) 

where fn^x is a parameter. The motivation of this ansatz fol- 
lows from the realization that if every atom in the universe 
is ionized exactly once, then in the beginning of reionization, 
while jCeff is sufficiently small, jCgff ~ N^/a-^ The factor fn^x 
therefore accounts for the loss of ionizing photons to recombi- 
nations and for more complicated dependence of Xeft on {x) at 
the later stages of reionization. This factor cannot be too small 
(or the universe would never be reionized), and any value for 
fn-,x above about 0.5 (half of photons lost for recombination) 
does not affect our conclusions. Therefore, in the rest of this 
paper we adopt f„^x = 0.75 as our fiducial value. In reality, 
f„^x must be a function of time, but since it is not likely to 
be much lower than 1, the exact time dependence of f„^x is 
unimportant at the level of precision of our approximations 
(which is dominated by the uncertainties on the escape frac- 
tion from galaxies and the lack of knowledge for the specific 
parameters of our quasar model). 

3. RESULTS 

In the rest of this paper we consider two different cosmolog- 
ical models. The first one (which we call "WMAP-5") is the 
most likely joint fit to the fifth year WMAP data, Baryonic 
Acoustic Oscillations from SDSS, and super novae, as pre- 
sented in the second column of Table 1 from iKomatsu et alJ 

* If we ignore the second term in equation |6] wliich is small for S> 
100 eV, and expand /si in Taylor series around x = 0, then for sufficiently 
small X, x^ff = {x) . 
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TABLE 1 
COSMOLOGICAL PARAMETERS 



Pai'ameter 


WMAP-5 


WMAP-5UP 




0.28 


0.28 
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0.82 


0.84 
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0.96 


0.98 




10-= r 
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Fig. 2. — The total ionizations-to-atom ratio as a function of redsliift for our 
fiducial model ("large seeds, high spins" in WMAP-5UP cosmology) for 3 
values of the mean energy of ionizing photons from quasars . The adopted 
reionization criterion, 1 < N^i„ < 3 at z = 6, is shown as a thick black segment 
with eiTor-bars. 



(I2009h . The second model ("WMAP-5UP") is obtained from 
the first one by increasing both the ampHtude of the density 
fluctuations erg and the scalar spectral index ns upward by 1 
standard deviation. The relevant values of cosmological pa- 
rameters for these two models are listed in Table[T] 

Figure |2] demonstrates the sensitivity of the total 
ionizations-to-atom ratio N^^ja on the assumed value for the 
mean energy of ionizing photons from quasars E^. Because 
ionizations from the quasar population are dominated by sec- 
ondary ionizations (which are independent of E^) for any 
plausible value of E^, this parameter has only a mild effect 
on our results. 

Figures |3]and|4]now present our main result: the computed 
N^la as a function of redshift for two cosmologies and various 
parameters of the quasar model. The criterion for the universe 
to be reionized at z = 6 (1 < N^i^ < 3) we adopt from Gnedi^ 
(j2008), where it is discussed and justified. 

Clearly, the WMAP5 model (with the most likely values of 
t78 and ns) is somewhat sh ort of the reion ization requirement. 
As has been discussed by iGnedinI (1200 8). this is not neces- 
sarily a serious problem, since the likely uncertainty on our 
extremely simple model is a factor of 2 or so. Never-the-less, 
slightly higher values for a% and/or ns provide a wider breath- 
ing space for reionization modeling. 

In order to illustrate a plausible estimate for the uncer- 
tainties jiue_to observational errors on the escape fraction 
(c.f. Shaplev et al. 2006), high redshift luminosity functions 
(IBouwens et al..,2007.) . etc, we show in Figure the "large 
seeds - high spin" quasar model in the WMAP-5UP cosmol- 
ogy (which we adopt as our "fiducial model") togeth er with 
the es timated errors. Notice, that formally Bouwe ns et al.l 
(|2007) galaxy luminosity function at z > 6 should be con- 
sidered as an upper limit; we emphasize this by having the 
hatched area in Fig.|5]unclosed from below. We also show the 
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Fig. 3. — The total ionizations-to-atom ratio as a function of redshift 
for different models for quasar contribution (red and blue lines) as well as 
contribution from the stars alone (black solid line) for the Eddington accre- 
tion model for two chosen cosmologies. The adopted reionization criterion, 
1 < N^i„ < 3 at z = 6, is shown as a thick black segment with error-bars. 



10 



0.1 



10- 



10 ■* 



I ' ' ' I ' ' ' I ' ' ' 
WMAP-5 Cosmology 

Merloni-Heinz Accretion 



Stars only 
Small seeds, low spins 
Small seeds, high spins 
Large seeds, low spins 
Large seeds, high spins 
I . I I ... I > I 



10 



12 



14 



10 



10- 



10- 



10- 



: 1 1 1 1 1 1 


1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 , 




WMAP-5UP Cosmology : 




Merloni-Heinz Accretion -= 


r stars 


only 


- Small 


seeds, low spins 


'_ Small 


seeds, high spins 


-. Large 


seeds, low spins E 


: Large 


seeds, high spins I 

, 1 , < , 1 < , , 1 , , , 1 , 



Fig. 4. — Same as Fig.|4] but for the Merloni & Heinz accretion model. 

evolution of the effective ionization fraction jCeff as an illustra- 
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z 

Fig. 5. — The total ionizations-to-atom ratio as a function of redsliift our 
fiducial model ("large seeds, high spins" in WMAP-5UP cosmology - red 
lines). The hatched area shows our estimate of the observational uncertainty 
(for a given theoretical model). The blue line shows the effective ioniza- 
tion fraction A'eff from equation (T)- Dotted, short-dashed, and long-dashed 
red lines show the contributions from primary and secondary ionization from 
quasars and from stars respectively. 

live (but highly approximate) reionization history. 

4. DISCUSSION 

Since our quasar models cover a wide range of physically 
plausible possibilities, we can draw some general conclusions 
from Fig. [3] and m In particular, we notice that (i) by z « 6 
the quasar contribution to the total ionization budget becomes 
mostly sub-dominant, not exceeding 50% in the best case, and 
likely falling below 20% or so. Never-the-less, at z « 8 in 
all our models quasars contribute from over 50% up to 90% 
of ionizations. This conclusion is particularly important to 
theoretical modeling and future observational measurements 
of the expected redshifted 21 cm emission from neutral hy- 
drogen in the reionization era. While most of the currently 



planned pathfinder experiments rely, in part, on the existing 
simulations for critical design decisions, the predictions from 
those simulations - none of which include a quasar contri- 
bution - become inaccurate at z > 8. It is, therefore, impor- 
tant to keep in mind that much theoretical work still needs 
to be done before even the basic observables for the incom- 
ing 21 cm experiments (like the fluctuation power spectrum) 
can be computed theoretically with ^20% precision (within a 
given cosmological model). 

Most of the contribution to the quasar ionizing budget is 
dominated by MBHs with mass < 10'' M©. Such small, low- 
luminosity MBHs do not contribute to the bright end of the 
luminosity function of quasars, and are therefore difficult to 
account from simple extrapolations of the luminosity func- 
tion of quasars. These small holes are not hosted in extremely 
massive galaxies residing in the highest density peaks (5 to 
6(7 peaks), but are instead found in more common, "normal" 
systems, ~ 3(t, peaks. Future generation of space-based tele- 
scopes, such as JWST and IXO, are likely to detect and con- 
strain the evolution of the population of accreting massive 
black holes at early times (z < 10). 
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